#ifndef _ADC_2_H
#define _ADC_2_H
#include<vector>
#include<eigen3/Eigen/Core>
#include"system_hf.h"

using Eigen::MatrixXd;
using Eigen::VectorXd;
using Eigen::Ref;

struct block
{
  int dim;
  int l,jj,tz;
};

struct state_2p1h
{
state_2p1h(int _p1,int _p2,int _h1,int _J):p1(_p1),p2(_p2),h1(_h1),J(_J){};
  int p1,p2;
  int h1;
  int J;///< p1 and p2 couple to J
};

struct state_2h1p
{
state_2h1p(int _h1,int _h2,int _p1,int _J):h1(_h1),h2(_h2),p1(_p1),J(_J){};
  int h1,h2;
  int p1;
  int J;//h1 and h2 couple to J
};

class ADC2
{
 public:
  typedef System_Tz_HF::OrbitalType HFOrbital;
  typedef vector<HFOrbital> HFOrbitals;
  typedef vector< vector<int> > OrbitalGroups;
  
 ADC2(const System_Tz_HF*Ptr):pSystem_Tz_HF(Ptr),HFOrbs(pSystem_Tz_HF->Orbitals),\
    Groups(pSystem_Tz_HF->Groups){}
  void setup();
  void setup2p1h();
  void setup2h1p();

  /// return the bra and ket onebody element in nb th group of orbitals.
  double get1B(int nb,int bra,int ket)
  {
    return pSystem_Tz_HF->get1B(pSystem_Tz_HF->Groups[nb][bra],pSystem_Tz_HF->Groups[nb][ket]);
  }
  ///(nb1,a1) represent a orbital
  double get2B(int nb1,int a1,int nb2,int a2,int nb3,int a3,int nb4,int a4,int J)
  {
    int a=pSystem_Tz_HF->Groups[nb1][a1];
    int b=pSystem_Tz_HF->Groups[nb2][a2];
    int c=pSystem_Tz_HF->Groups[nb3][a3];
    int d=pSystem_Tz_HF->Groups[nb4][a4];
    return pSystem_Tz_HF->get2B(a,b,c,d,J);
  }
  
  void solveDysonMatrix();

  void solveDysonMatrix_Lanczos(int p,int q);
  void buildLanczosMatrix(Ref<MatrixXd> A,Ref<MatrixXd> C,const MatrixXd & M_A,const MatrixXd & M_C,const vector<int>& n_iter);
  ///one self-consistent iteration
  void solveDysonMatrix_Lanczos_iter(int p,int q);
  void buildLanczosMatrix_iter(Ref<MatrixXd> A,Ref<MatrixXd> C,const MatrixXd & M_A,const VectorXd & M_C,const vector<int>& n_iter);
  
  void calEnergy();

  void printResults() const;
  
  double sep_energy;
  double Energy;
  double num;
  
  //the outside vector for different block
  vector<block> blocks;
  vector< vector<state_2p1h> > states_2p1h;
  vector< vector<state_2h1p> > states_2h1p;
  vector< MatrixXd > Z;
  vector< VectorXd > poles;
  const System_Tz_HF* pSystem_Tz_HF;
 private:
  const HFOrbitals & HFOrbs;
  const OrbitalGroups & Groups;
  int totalOrbitals;
  int FermiSurface;
  int numGroups;
};

#endif

